Hazard assessment and hydrodynamic, morphodynamic, and hydrological response to Hurricane Gamma and Hurricane Delta on the northern Yucatán Peninsula

Barrier islands in tropical regions are prone to coastal flooding and erosion during hurricane events. The Yucatán coast, characterized by karstic geology and the presence of barrier islands, was impacted by Hurricane Gamma and Hurricane Delta in October 2020. Inner shelf, coastal, and inland observations were acquired simultaneously near a coastal community (Sisal, Yucatán) located within 150 km of the hurricanes' tracks. In the study area, Gamma moved slowly and induced heavy rain, mixing in the shelf sea, and strong winds (>20 m s⁻¹). Similar wind and wave conditions were observed during the passage of Hurricane Delta; however, a higher storm surge was measured due to wind setup and the drop (<1000 mbar) in atmospheric pressure. Beach morphology changes, based on GPS measurements conducted before and after the passage of the storms, show alongshore gradients ascribed to the presence of coastal structures and macrophyte wracks on the beach face. Urban flooding occurred mainly on the back barrier associated with heavy inland rain and the coastal aquifer's confinement, preventing rapid infiltration. Two different modeling systems, aimed at providing coastal flooding early warning and coastal hazard assessment, presented difficulties in forecasting the coastal hydrodynamic response during these seaward-traveling events, regardless of the grid resolution, which might be ascribed to a lack of terrestrial processes and uncertainties in the bathymetry and boundary conditions. Compound flooding plays an important role in this region and must be incorporated in future modeling efforts.

1 Introduction

Barrier islands are highly vulnerable to terrestrial, atmospheric, and oceanic drivers associated with hurricane events and global climate change effects (Irish et al., 2010; Zinnert et al., 2017). The mean sea level rise rate is accelerating, and the proportion of major tropical cyclones has increased over recent decades due to the effects of climate change (Knutson et al., 2020). Therefore, knowledge of coastal dynamics during hurricane events and their impact on barrier islands is important to understand such coastal ecosystems' natural resistance and resilience.

Significant advances have been achieved in recent decades regarding hurricane research. While Emanuel (2021) found an increased frequency of tropical cyclones in the North Atlantic, there is no clear trend in the increase in tropical cyclone frequency due to climate change. Nevertheless, most studies find an increased proportion of the most extreme events (i.e., categories 4 and 5 on the Saffir–Simpson scale) in the context of climate change (Knutson et al., 2020), increasing associated hazards by the second half of the century. Furthermore, a poleward migration of the location of the maximum lifetime intensity of tropical cyclones has been found (Kossin et al., 2014), as well as an increase in rapid intensification (Bhatia et al., 2019; Emanuel, 2017), increasing the hazards from tropical cyclones in higher latitudes and hence representing a challenge for emergency management.

Field observations provide valuable information to improve our understanding of the drivers and coastal response during extreme events (e.g., Valle-Levinson et al., 2002; Du et al., 2019; Mieras et al., 2021). Previous studies have demonstrated the important role of compound flooding driven by rain, storm surge, and groundwater processes (e.g., Wahl et al., 2015; Valle-Levinson et al., 2020; Housego et al., 2021). These studies are required to improve coastal hazard modeling and implement mitigation measures at both local and regional scales. However, concurrent measurements of meteorological, oceanic, coastal, and hydrological processes during extreme events are scarce in tropical systems; hence further research is warranted. Moreover, recent studies have pointed out the need to investigate storm impacts from an interdisciplinary point of view. Considering the importance of numerical modeling for developing early warning systems and emergency management plans, an improvement in the understanding of the feedback between terrestrial–atmospheric–oceanic processes needs to be incorporated into numerical models to estimate their impacts in coastal areas. This is a challenging task for a low-lying (karstic) coast where the groundwater aquifer discharges toward the shore.

The Yucatán Peninsula is located in an area of high tropical cyclone activity. However, due to its geographic location, the northern coast (facing the Gulf of Mexico) is less prone to direct hurricane landfalls than other regions in the Gulf of Mexico and the western Caribbean Sea. For instance, of the 163 tropical storm events from 1842 to 2020, 64 % landed on the eastern coast facing the Caribbean Sea (Rivera-Monroy et al., 2020), compared to 35 % on the northern coast. Human settlements on a barrier island along the northern Yucatán coast are mainly devoted to artisanal fisheries (Paré and Fraga, 1994) and recreational beach houses (Meyer-Arendt, 2001). These coastal communities are highly vulnerable to storm events due to their high exposure and potential communication breakdown with the inland due to high water levels.

Concurrent observations of atmospheric, terrestrial, and oceanic effects and coastal impacts during tropical cyclones are rare or inexistent in the region. Thus, this study investigates the effects and impacts of hurricane passages, specifically the tropical cyclones Delta and Gamma, on a barrier island located on the northwestern Yucatán coast for the first time. The observations of this coastal community are relevant as they represent the environmental, morphological, anthropogenic, and ecological conditions prevailing in this region (i.e., northern Yucatán Peninsula).

The outline of this paper is the following. An overview of the study area and the meteorological events (Gamma and Delta) is provided in Sect. 2. Section 3 describes the materials and methods employed in this work, including the data acquisition and analysis, and the implementation of the numerical models. Section 4 presents the observations of atmospheric, oceanic, and hydrological conditions associated with tropical storms and their impact on the coast of Sisal (Yucatán). Moreover, the capabilities and limitations of hydrodynamic numerical models for simulating waves and extreme water levels are investigated. Finally, discussions (Sect. 5) and concluding remarks (Sect. 6) are presented.

2 Study area

The study area is the town of Sisal, Yucatán, located on the southeastern Gulf of Mexico coast (21^∘09^′56.20^′′ N, 90^∘02^′26.44^′′ W). Sisal is a small community situated on a barrier island on the northwestern Yucatán Peninsula (Fig. 1a), 50 km from the city of Mérida, with a population of less than 2000, dedicated mainly to fishing activities and, more recently, to eco-tourism. The main infrastructure found in this community is a sheltered port devoted to artisanal fisheries and a local road passing through the wetland, connecting the barrier island with the hinterland. The area is of high ecological importance due to its biological diversity and because it is surrounded by two natural parks, Ciénegas y Manglares de la Costa Norte de Yucatán and the state reserve of El Palmar to the west. Emblematic species such as jaguars (Panthera onca), crocodiles (Crocodylus moreletii), flamingos (Phoenicopterus ruber), and sea turtles (Eretmochelys imbricata and Chelonia mydas) are found in the tropical dry forest, wetlands, and beaches surrounding this coastal town.

Figure 1 Study area showing (a) best track positions for Hurricane Gamma and Hurricane Delta (National Hurricane Center, https://www.nhc.noaa.gov/data/tcr/index.php?season=2020&basin=atl, last access: December 2022), (b) the location of the acoustic Doppler current profiler (ADCP) and the monitoring wells, and (c) coastal monitoring systems and beach transects.

The study area is characterized by intense sea breeze events (Figueroa-Espinoza et al., 2014) and winter storms (fall–winter) associated with Central American cold surge events (Medina-Gómez and Herrera-Silveira, 2009; Kurzcyn et al., 2021). Typical winter storm wave conditions reach a significant wave height (H_s) >2 m and a peak wave period (Tp) >8 s from the NNW. Sea breezes drive low-energy high-angle short-period waves (H_s<1 m, Tp=3.5 s, NE) that drive a persistent westward alongshore current (Torres-Freyermuth et al., 2017). Across the Yucatán shelf, winds and mesoscale circulation drive the currents mainly toward the west (Enriquez et al., 2010; Torres-Freyermuth et al., 2017). The mean sea level in this area presents a seasonal variability ascribed to along-shelf currents on the western Gulf of Mexico and low-frequency atmospheric pressure variability (Zavala-Hidalgo et al., 2003). The tidal regime in the study area is diurnal micro-tidal, with a spring tidal range of 0.75 m (Valle-Levinson et al., 2011).

The 2020 hurricane season was the most active season on record in the North Atlantic basin (Blunden and Boyer, 2021). Hurricane Gamma formed on 2 October in the western Caribbean Sea as a tropical depression southeast of Cozumel, Mexico. The hurricane made landfall on the eastern Yucatán Peninsula coast near Tulum, Mexico, on 3 October and then weakened into a tropical storm while crossing the peninsula and reaching the Gulf of Mexico via the northern coast (Fig. 1). Gamma interacted with the circulation associated with the formation of Hurricane Delta and moved southwestward to make landfall near Nichili, Mexico, and dissipated on 6 October (Latto, 2021). Hurricane Delta formed from a tropical wave in the Atlantic and attained the category of major hurricane on 6 October, before undergoing rapid intensification and weakening before landing. It made landfall on 7 October, on the northeastern portion of the Yucatán Peninsula near Puerto Morelos, Mexico, around 10:30 UTC (20.848^∘ N, 86.875^∘ W). It continued its path inland and moved to the southern Gulf of Mexico by 18:00 UTC on 7 October with winds of 38 m s⁻¹, reaching the coast of Dzilam de Bravo (Yucatán) (21.393^∘ N, 88.892^∘ W) as a category-1 hurricane (160 km h⁻¹). During its pass across the northern Yucatán Peninsula, it dumped 50–100 mm of rain. Delta made landfall near Creole, Louisiana, on 9 October as a category-2 hurricane and weakened to become an extratropical cyclone (Cangialosi and Berg, 2021).

3 Materials and methods

3.1 Data acquisition

Monitoring systems installed in the study area were employed to characterize the atmosphere, the ocean, beach morphology, and the coastal aquifer response to the storm forcings (Fig. 1b and c and Table 1) west of the hurricane tracks.

Table 1 Instruments and measured variables.

A meteorological station, located at 10 m height and 100 m from the shoreline, measured the wind intensity and direction, the air temperature, and the atmospheric pressure at 10 Hz. A tidal gauge inside the port of Sisal recorded the mean sea level every 1 min with an ultrasonic sensor. An acoustic Doppler current profiler (ADCP), deployed 10 km offshore at 11 m depth, measured the wave parameters (H_s; Tp; θ: mean wave angle), the current profile, the sea surface height (η), and the sea bottom temperature (T_b). Waves were measured by taking 2048 samples at 2 Hz every hour, while the rest of the ADCP data were obtained every 20 min, averaging the first 60 s of the observations.

A beach monitoring program has been conducted regularly to investigate beach morphodynamics in the study area (Medellín and Torres-Freyermuth, 2019, 2021; Franklin et al., 2021). For this study, pre-storm (30 September 2020) and post-storm (14 October 2020) beach surveys were conducted along 40 equally spaced cross-shore transects located east (updrift, P01–P20) and west (downdrift, P21–P40) of the Sisal port (see Fig. 1c), encompassing a 4 km stretch of coast. Differential global positioning systems (DGPSs) with real-time kinematics (RTK) were employed to conduct the beach survey. A reference station is located at a fixed location (top of a building), and the rover is carried on a backpack. Ground control points were measured at the beginning and end of each survey to correct the rover height. The DGPS RTK measurements have a horizontal and vertical accuracy of 0.010 and 0.020 m, respectively, as reported by the manufacturer. Additionally, a test was carried out to assess the method's accuracy by comparing elevation measurements taken along a straight line on a parking lot with the GPS rover fixed on a pole with a bubble level against carrying the GPS rover in a backpack by two different users. Maximum differences in the vertical measurements were less than 0.04 m. Therefore, we expect that the maximum error associated with estimating the beach profiles' elevation is of that magnitude. The beach profiles started behind the foredune and reached a water depth of approximately 0.5–1.5 m, depending on existing wave conditions, the tidal level, and the presence of macrophytes.

A video monitoring system, placed at a height of 43 m, located 300 m west of the jetty and 100 m inland, acquired time exposure images over 10 min at 7.5 Hz of 2 km along the coast every 30 min during daytime hours (Arriaga et al., 2022). A gap in the images occurred between 6 and 19 October due to power failure. The spatial resolution of the different components of this camera system has been previously described by Arriaga et al. (2022). Given that here we are interested mainly in coverage, this result has to be translated into areas. For example, along the shoreline near the cameras (300 m easting), a pixel translates to an area of 0.01 m², whereas at 1400 m easting a pixel covers an approximate size of 2.5 m². Finally, at the farthest point, near the pier, the resolution is on the order of 10 m². The calibration performed following Simarro et al. (2017), particularly the simplified mathematical model referred to as M2 in Simarro et al. (2020), resulted in an average reprojection error of 0.7 pixels.

The coastal aquifer response was characterized by groundwater pressure, temperature, and salinity measured every 30 min using pressure transducers (HOBO) installed in three monitoring wells (Fig. 1b). The well W7a is located close to the Sisal port, whereas W5 and W4 are located 5 and 20 km inland, respectively. A detailed description of the monitoring wells can be found in Canul-Macario et al. (2020).

3.2 Data analysis

Time series analyses, obtained from different in situ instruments, were carried out to characterize the main drivers and responses.

3.2.1 Monitoring wells

Pressure measured at the monitoring wells is converted to hydraulic head as follows: atmospheric pressure, p_atm, is subtracted from the pressure, p, measured at the well (kPa) and then converted to water column height H (m), using a conversion factor of 0.102 using Eq. (1). Finally, the measurements are referenced to the same datum using the elevation of a known point at the well casing, z_well, and one measurement from this point to the water depth, H₀, at a known time, t₀. The equations used are the following:

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }H\left(t\right)=\mathrm{0.102}\left(p\left(t\right)-p_{\mathrm{atm}}\left(t\right)\right),\text{(2)}& {\displaystyle }& {\displaystyle }h\left(t\right)=H\left(t\right)-H\left(t_{\mathrm{0}}\right)+\left(z_{\mathrm{well}}-H_{\mathrm{0}}\right),\end{array}$$

where h(t) is the hydraulic head (m) time series that represents the water level referenced to a given datum for unconfined aquifers. To obtain the relative water head, the mean value is subtracted from the time series.

3.2.2 Acoustic Doppler current profiler

Bulk wave statistics (H_s, Tp, θ) were computed from the 2048 s burst intervals using the PUV method. For the ADCP measurements, the tidal signature was removed from the current profile and the sea surface height by applying a low-pass Lanczos filter, eliminating frequencies ≥ $\mathrm{1}/\mathrm{48}$ h. Ocean currents were then referenced to the angle of maximum variance, orienting them 25^∘ counterclockwise from the east. To compare the atmospheric and oceanic measurements, these were homogenized in time by taking the daily averages of each observation.

Heat fluxes were computed as follows. Sensible heat (Q_h in W m⁻²) was based on Gill (1982) and Talley et al. (2011):

$$\begin{array}{}\text{(3)}& Q_{\mathrm{h}}={\mathit{\rho }}_{\mathrm{a}}{C}_{\mathrm{p}}{C}_{\mathrm{h}}\mathrm{wVel}\left[\mathrm{SST}-\left(T_{\mathrm{air}}+\left({\displaystyle \frac{\mathrm{9.8}}{\mathrm{1000}}}\right)z_{\mathrm{air}}\right)\right],\end{array}$$

where ρ_a is the air density (in kg m⁻³), C_p is the specific heat capacity of air at constant pressure (in J kg⁻¹ K⁻¹), C_h is the bulk sensible heat transfer coefficient, wVel is the wind velocity (in m s⁻¹), SST is the sea surface temperature (from satellite remote sensing in ^∘C), T_air is the air temperature, and z_air is the height where the air temperature was taken (∼10 m a.s.l. – above sea level). C_h was made dependent on the wind velocity and $\mathrm{\Delta }\mathrm{SST}=\mathrm{SST}-T_{\mathrm{air}}$ according to table values by Smith (1988).

Latent heat (Q_e in W m⁻²) was estimated following Gill (1982) and Castro et al. (1994):

$$\begin{array}{}\text{(4)}& Q_{\mathrm{e}}={\mathit{\rho }}_{\mathrm{a}}{C}_{\mathrm{e}}\mathrm{wVel}{L}_{\mathrm{v}}\left(q_{\mathrm{s}}-q_{\mathrm{a}}\right),\end{array}$$

where C_e is the exchange coefficient, wVel is the wind velocity (in m s⁻¹), q_s is the saturation specific humidity of the sea surface, q_a is the specific humidity of air, and L_v is the latent heat of evaporation (in W m⁻²). C_e was made dependent on the wind velocity and ΔT according to table values by Bunker (1976), whereas the specific humidity saturation of the sea surface (q_s) was calculated with (Gill, 1982)

$$\begin{array}{}\text{(5)}& q_{\mathrm{s}}={\displaystyle \frac{\mathrm{0.62197}{e}_{\mathrm{w}}}{\left(P_{\mathrm{atm}}-\mathrm{0.378}{e}_{\mathrm{w}}\right)}},\end{array}$$

where e_w and q_a are explained above. The latent heat of evaporation (L_v) was measured as

$$\begin{array}{}\text{(6)}& L_{\mathrm{v}}=\mathrm{2.5008}\times {\mathrm{10}}^{\mathrm{3}}-\mathrm{2.3}\mathrm{SST}.\end{array}$$

Moreover, from the weather station data, the impact of the hurricane winds on the water column was analyzed by computing the Ekman surface velocity (U_E) and the Ekman layer depth (D_e). U_E was estimated following Rio et al. (2014):

$$\begin{array}{}\text{(7)}& U_{\mathrm{E}}\left(z\right)=\mathit{\beta }\left(z\right)\mathit{\tau }{e}^{\left(i\mathit{\theta }\right(z\left)\right)},\end{array}$$

where the parameters β and θ at the surface are β(0)=0.61 and θ(0)=30.75, respectively. τ is the surface wind stress (in N m⁻²):

$$\begin{array}{}\text{(8)}& \mathit{\tau }=C_{\mathrm{d}}{\mathit{\rho }}_{\mathrm{a}}{\mathrm{wVel}}^{\mathrm{2}}.\end{array}$$

C_d is the drag coefficient; wVel is the weather station wind velocity (in m s⁻¹); and ρ_a is the air density defined above. The Ekman layer depth (in m) was based on Cushman-Roisin and Beckers (2011):

$$\begin{array}{}\text{(9)}& D_{\mathrm{E}}=\mathrm{0.4}{\displaystyle \frac{u^{*}}{f}},\end{array}$$

where f is the Coriolis parameter estimated at the position of the ADCP and u^* is the turbulent velocity:

$$\begin{array}{}\text{(10)}& u^{*}=k{\displaystyle \frac{\mathrm{cVel}}{\log \left(\frac{z}{z_{\mathrm{0}}}\right)}}.\end{array}$$

Here, cVel (in cm s⁻¹) is the vertically averaged current velocity, k is the von Kármán constant (0.41), z is the mean current measurement depth (5.5 m in our case), and z₀ is the size of the ripples or gravel on the seafloor (∼0.05 m).

3.2.3 Satellite imagery

For the analysis of (1) surface winds, (2) sea surface salinity (SSS), and (3) sea surface temperature (SST), the following remote sensing images were used:

• i.

  MetOp Advanced SCATterometer (ASCAT) $\mathrm{1}/\mathrm{4}$^∘ daily wind and wind stress maps were used from the Centre d'Exploitation et de Recherche SATellitaire (CERSAT), at IFREMER, Plouzané (France). More details on the data, objectives, method, and computation algorithm are found in Bentamy and Croizé-Fillion (2012). Data and documentation are freely distributed at ftp://ftp.ifremer.fr/ifremer/cersat/products/gridded/MWF/L3/ASCAT/Daily/ (last access: January 2022).

• ii.

  The L3_DEBIAS_LOCEAN_v5 $\mathrm{1}/\mathrm{4}$^∘ 4 d sea surface salinity maps (Boutin et al., 2018) have been produced by the LOCEAN IPSL (UMR CNRS–UPMC–IRD–MNHN) laboratory and the ACRI-ST company that participate in the Ocean Salinity Expertise Center (CEC-OS) of the Centre Aval de Traitement des Données SMOS (CATDS). This product is distributed by the Ocean Salinity Expertise Center (CEC-OS) of the CNES IFREMER Centre Aval de Traitement des Données SMOS (CATDS), at IFREMER, Plouzane (France) (ftp://ext-catds-cecos-locean:catds2010@ftp.ifremer.fr/, last access: January 2022).

• iii.

  The NOAA $\mathrm{1}/\mathrm{4}$^∘ Daily Optimum Interpolation Sea Surface Temperature (Reynolds et al., 2002) was provided by the NOAA/OAR/ESRL Physical Sciences Division (PSD), Boulder, Colorado, USA, at https://www.esrl.noaa.gov/psd/ (last access: January 2022).

Surface winds, SSS, and SST observations were spatially interpolated on a 25 km radius from the ADCP location; salinity and temperature data were later transformed to conservative temperature and absolute salinity to estimate σ_θ (in kg m⁻³), based on the thermodynamic equation of seawater TEOS-10 (McDougall and Barker, 2011). The precipitation brought by these storms caused changes in the seawater density (σ_θ in kg m⁻³), which were inspected employing the SSS and SST data, whereas the remotely sensed surface wind stress was used to estimate the Ekman pumping (W_E in m s⁻¹), following the formulas proposed by Cushman-Roisin and Beckers (2011):

$$\begin{array}{}\text{(11)}& W_{\mathrm{E}}={\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}\left[{\displaystyle \frac{\mathrm{d}}{\mathrm{d}x}}\left({\displaystyle \frac{{\mathit{\tau }}^y}{f}}\right)-{\displaystyle \frac{\mathrm{d}}{\mathrm{d}x}}\left({\displaystyle \frac{{\mathit{\tau }}^x}{f}}\right)\right],\end{array}$$

where W_E (in m s⁻¹) is the vertical velocity estimated from the wind stress components (τ^x, τ^y) of the satellite surface winds near the ADCP location, f is the Coriolis parameter, and ρ₀ (in kg m⁻³) is the water density. Multiplying the Ekman pumping by the Ekman layer, the vertical advective transport (W_aT in m² s⁻¹) was obtained:

$$\begin{array}{}\text{(12)}& W_{\mathrm{aT}}=W_{\mathrm{E}}{D}_{\mathrm{E}}.\end{array}$$

3.2.4 Beach surveys

The most notable coastal impacts associated with the passage of storms are beach erosion and flooding. Beach profiles were employed to estimate the subaerial beach volume change by integration of beach elevation (z≥0 m) with respect to the cross-shore distance. The subaerial beach volume change can be readily obtained by subtracting the pre-storm from the post-storm beach survey. On the other hand, shoreline position was obtained by tracking the cross-shore location corresponding to z=0 m for each transect; hence shoreline change was estimated as the difference between the pre- and post-storm shoreline location. To estimate the coastal flooding during the peak of the storm, beach elevation changes between subsequent surveys were employed as a proxy for the maximum water levels at each transect.

3.2.5 Video camera system

Time exposure images were employed to observe the storm effects on macrophyte wrecking and dune vegetation and to estimate post-storm inundated areas. To quantify the impact in meters, image pixel coordinates were transformed to UTM coordinates by relating ground control points (GCPs) to pixel positions (Simarro et al., 2017). Unfortunately, the strong winds moved the orientation of the cameras and even the GCPs. To solve this, a previous image with known GCPs was used as a reference to stabilize the images of interest (Arriaga et al., 2022). Following the methodology of Rutten et al. (2021) to detect Sargassum on images, the support vector machine (SVM) algorithm was employed to detect the coverage of wrack, vegetation, and flooded areas (1, 2, 5, 20 October).

3.3 Numerical modeling

Two different numerical approaches were implemented in the study area. The numerical models were forced with different wind information and grid resolutions. The first approach aimed to forecast the atmospheric and oceanic conditions generated during the passing of the events, while the second approach aimed to determine the wave and storm surge hazards created by the resulting waves and storm surge. A description of each modeling approach and its implementation is provided below.

3.3.1 Forecast modeling

The numerical models WRF, WWIII, and ADCIRC were employed to forecast nearshore hydrodynamics as described below.

WRF model

The Weather Research and Forecasting (WRF V.3.9) model, developed by the National Center for Atmospheric Research (NCAR), is characterized by being compressible and non-hydrostatic, with terrain-following hydrostatic pressure vertical coordinates and Arakawa C horizontal grid staggering (Arakawa and Lamb, 1977). The model used the Runge–Kutta second- and third-order time integration schemes and the second- to sixth-order advection schemes in both the horizontal and the vertical. Moreover, it also uses a small time-split step for acoustic- and gravity-wave modes. For further details, refer to Skamarock et al. (2008). The operational forecasting system established at the Institute of Atmospheric Sciences and Climate Change (ICAyCC) at the National Autonomous University of Mexico (UNAM; Ocean-Atmosphere Interaction Group, 2020) was used. The physics model parameterizations are the Kain–Fritsch cumulus parameterization scheme (Kain, 2004), the RRTM (rapid radiative transfer model) scheme for longwave radiation (Mlawer et al., 1997), the Dudhia scheme for shortwave radiation (Dudhia, 1989), and the Yonsei University (YSU) scheme for the boundary layer (Skamarock et al., 2008; Hong et al., 2006). In addition, the five-layer thermal diffusion land surface model (LSM) was used. This scheme, although simple, is adequate for most mesoscale studies and estimates the energy balance at a low computational cost (Dudhia, 1996). The land use and land cover soil category map data used were obtained from the USGS database with 24 classes (Loveland et al., 2000). Here, the forecast employed two one-way nested computational domains. The first domain (D01) has a 15 km horizontal grid resolution and includes Mexico, the Gulf of Mexico, part of the Caribbean Sea, and part of the central Pacific. The second domain (D02) has a resolution of 5 km and includes the central part of the Mexican territory. The forecast employed 30 vertical levels in a log-normal distribution, with the top of the atmosphere fixed at 50 mbar. The model equations were integrated every 120 s. For the initial and boundary conditions, the numerical model was initialized with the Global Forecast System (GFS) model at 00:00 UTC data, every 6 h with a 1^∘ spatial resolution. The operational system produces a 5 d forecast; however for this work only the 24, 48, and 72 h forecasts were considered. It is known that a tropical cyclone's track and intensity forecast degrades rapidly, so it is not reliable beyond 72 h. Despite the performance of the forecasts being acceptable for this particular case (see Table 1), it was not considered for discussion. The correlation coefficients halved for wind, waves, and sea level parameters beyond 72 h.

WAVEWATCH III

The WAVEWATCH III (WWIII V.5.16) model is a third-generation model developed by NOAA (National Centers for Environmental Prediction, NCEP), characterized for solving the random-phase spectral action density balance equation for wave number direction spectra. The implicit assumption of this equation is that properties of the medium (water depth and current) and the wave field itself vary over timescales and space scales that are much larger than the variation scales of a single wave. This model also considers options for extremely shallow water (surf zone), as well as wetting and drying of grid points. The total wave energy and the local and instantaneous spectrum of the waves can be obtained as model outputs, where the latter can be reduced to a two-dimensional function. The parameters of significant wave height, mean period, and direction of propagation are also model outputs. The parameterizations of the operational forecasting system used in this work are the Cavaleri and Rizzoli (1981) term to represent the linear growth of the waves and the parameterization described by Tolman and Chalikov (1996) in the terms that define the integral growth of waves. To represent nonlinear processes, the discrete iteration approximation described by Hasselmann et al. (1985) was used. The bottom friction was represented with an empirical linear function described in Hasselmann et al. (1973). The forecast uses a rectangular, rectilinear grid. The model was implemented in two one-way nested computational domains: the world ocean, on one side the Pacific Ocean and, on the other, the Gulf of Mexico. Atmospheric forcing was obtained from the Global Forecast System (GFS from NCEP; https://www.ncdc.noaa.gov/data-access/model-data/model-datasets/global-forcast-system-gfs, last access: October 2020) and the Weather Research and Forecasting (WRF) model, respectively, for each domain. The bathymetry used was the ETOPO1 elevation database from the National Centers for Environmental Information (NOAA; https://www.ngdc.noaa.gov/mgg/global/, last access: January 2020), with a spatial resolution of 1 arcmin. The operational system produces a 5 d, 3-hourly forecast; however, only the 24, 48, and 72 h forecasts were considered for this work.

Table 2 Setup characteristics for the WRF, WWIII, and ADCIRC models.

ADCIRC

The ADvanced CIRCulation Model for Oceanic, Coastal and Estuarine Waters (ADCIRC V.52.30.13 with NetCDF file support) is a system of programs for solving time-dependent free surface circulation and transport problems in 2D and 3D. These programs solve the movement equations for a rotating fluid through the Boussinesq and hydrostatic pressure approximations, both discretized in space by the finite-element method and in time by the finite-difference method. This way of solving the movement equations allows the use of highly flexible unstructured grids. ADCIRC calculates the surface elevation from the generalized wave-continuity equation (GWCE) and the current velocity from the momentum equations. All nonlinear terms have been retained in these equations. ADCIRC applications include wind and tidal circulation modeling, flood and storm surge analysis, dredging and disposal feasibility studies, and larval transport studies, as well as for nearshore marine operations. The configuration used in this work has a resolution equal to or less than 500 m along the Mexican coast, reducing the resolution to 4 km for the northern Gulf of Mexico. Along the open boundary, eight harmonic tidal constituents are prescribed (M2, S2, K2, N2, K1, O1, P1, and Q1) obtained from TPXO (https://www.tpxo.net/global/tpxo9-atlas, last access: January 2020; Egbert and Erofeeva, 2002). Surface forcings, including hourly winds and sea level atmospheric pressure, were taken from the WRF atmospheric model described in the previous section. Initial conditions were obtained from the atmospheric model output after regridding. Table 2 summarizes the setup for each of the aforementioned models.

3.3.2 Hazard assessment modeling

The numerical models MIKE 21 HD FM (hydrodynamic) and MIKE 21 SW (waves) were employed to assess wave and storm surge conditions. MIKE 21 SW is a third-generation spectral wave model which resolves the wave action equation employing non-structured grids using a finite volume. MIKE 21 HD FM (DHI, 2022) solves the Reynolds-averaged Navier–Stokes (RANS) equations with the Boussinesq and hydrostatic pressure approximations. The numerical model employs flexible meshes, allowing the resolution in the area of interest to be increased using a finite volume. This numerical model assumes a Coriolis force, baroclinic density, eddy viscosity of 0.28 based on the Smagorinsky formulation, and constant wind friction coefficient (0.001255 for wVel < 7 m s⁻¹ and 0.002425 for wVel > 7 m s⁻¹). MIKE 21 SW is implemented in stationary mode with a logarithmic discretization in the frequency domain and directional spectra divided into 32 bins. The wave and hydrodynamic models employ a computational domain covering the Gulf of Mexico with an 8 km resolution near the coast. Moreover, the resolution in the hydrodynamic model further increases in coastal areas up to 40 m in the area of interest. The topographic data used were obtained from a 1 m resolution lidar survey from 2011 of the Yucatán coast, while the bathymetry was obtained from local surveys complemented with ETOPO1 data (Amante and Eakins, 2009).

4 Results

4.1 Hurricane Gamma and Hurricane Delta

4.1.1 Atmospheric conditions

The passage of Hurricane Gamma induced an increase in the wind speed at Sisal, Yucatán, reaching a peak magnitude of 20 m s⁻¹ on 3 October 2020 and maintaining sustained winds around 15 m s⁻¹ for the following days until 6 October (Fig. 2a). The wind direction switched from the NE to the NNW on 4 to 6 October. The wind velocity dropped to 5 m s⁻¹ on 6 October, increasing suddenly to 20 m s⁻¹ on 7 October due to the passage of Delta. A sustained drop followed the wind increase throughout the same day due to Delta's fast translation speed. The wind direction switched from the NW to the SW as the storm passed. The atmospheric pressure (Fig. 2b) shows a significant decrease below 1000 mbar during the passage of Delta. The atmospheric temperature remained relatively steady, around 25 ^∘C, and the diurnal variability associated with the sea breeze recovered after 8 October (Fig. 2c). Heavy rain (50 mm) fell during 2 October followed by a lower peak on 7 October (Fig. 2d).

Figure 2 (a) Wind magnitude and direction, (b) atmospheric pressure, (c) air temperature, and (d) precipitation measured by the weather station at Sisal, Yucatán.

4.1.2 Oceanic conditions

In situ measurements at 11 m water depth allowed us to assess the effects of hurricanes Gamma and Delta on the nearshore sea. The ADCP time series spanned more than 3 years of data. However, Fig. 3 focuses on the daily averages of different variables from July to October 2020 to highlight the effects of hurricanes Gamma (5 October) and Delta (7 October). Figure 3a shows the wind velocity and the sea surface height oscillation before, during, and after the passage of such atmospheric events. The winds accelerated to reach 14.4 m s⁻¹ (52 km h⁻¹), blowing from the north and promoting a sea level setup (Δη) of 9 cm associated with Gamma and a 30 cm setup during Delta. Ocean currents increased to reach the maximum value that has been registered in the time series since 2018 (59 cm s⁻¹) during Gamma (Fig. 3b).

Figure 3 Time series of (a) wind velocity and sea surface height (in blue), (b) current velocity and significant wave height (in red), (c) Ekman currents and Ekman layer depth (in red), (d) sea bottom temperature and air temperature (in red), (e) seawater density and precipitation (in red), and (f) sensible and latent heat (in red), estimated from the meteorological station, the moored ADCP, and the satellite data. The dates of the passage of storms Gamma and Delta are highlighted. The date format is dd/mm/yyyy.

Figure 4 Mode 1 of the empirical orthogonal function (EOF) analysis for the ADCP alongshore currents: (a.1) temporal distribution and (a.2) spatial distribution. (b) Time–depth distribution of the alongshore current. The dates of the passage of storms Gamma and Delta are highlighted (dashed black and solid orange lines). Blue (red) colors represent a westward (eastward) flow. The date format is dd/mm/yy.

Furthermore, we analyzed the influence of wind stress over the nearshore sea by looking at the Ekman layer depth (D_E) and the surface currents (Ekman currents, U_E). The D_E average value for the ADCP location is 70±46 m (estimated from time series since 2018); that is, the Ekman layer commonly encompasses the entire water column (Table 3). In general, regional surface winds can generate an Ekman layer depth greater than 10 m 92 % of the time. During the passage of Delta and Gamma, the Ekman layer depth attained a mean value of 169 m and reached a maximum value of 350 m, implying significant water mixing throughout the water column and sediment resuspension due to seabed friction. U_E accelerated to the maximum value (26 cm s⁻¹) during Gamma. A comparison between the maximum current magnitude and the maximum currents due to the surface wind stress (U_E) suggests that the latter represented 44 % of the current magnitude (Fig. 3c). The vertical advective transport (W_aT) generated by Gamma and Delta (time series not shown) is limited by the depth of the study region, which restricts the Ekman layer depth. The vertical advective transport showed an average positive value during Gamma ($\mathrm{1.6}\times {\mathrm{10}}^{-\mathrm{4}}$ m² s⁻¹), i.e., an upwelling transport, and a sudden change during Delta from positive to negative values, changing from upwelling to downwelling transport in 2 d (from $\mathrm{4.0}\times {\mathrm{10}}^{-\mathrm{4}}$ to $-\mathrm{6.2}\times {\mathrm{10}}^{-\mathrm{4}}$ m² s⁻¹). At the same time, Gamma and Delta promoted a drop in air temperature (T_air) and sea bottom temperature (T_b) of 5 and 3 ^∘C, respectively (Fig. 3d). Moreover, the accumulation of freshwater due to the high precipitation volumes induced a decrease in the sea surface density, exhibiting low values (1022.7 kg m⁻³) at the end of these storms (Fig. 3e).

Table 3 Daily mean and maximum values during Gamma and Delta.

To investigate the heat exchange between the sea surface and the atmosphere, we estimated the sensible and the latent heat fluxes (Fig. 3f). The sensible heat (Q_h) is a proxy for the heat gain or loss from the sea surface due to thermal gradients between this and the adjacent air. On the other hand, the latent heat (Q_e) represents the heat exchange ascribed to evaporation/condensation processes between the sea surface and the atmosphere. A positive (negative) value represents a heat output (input) from the sea. Figure 3f shows the maximum sensible heat loss during the passage of Gamma. Thus, field observations suggest that both Gamma and Delta absorbed heat from the sea, although Delta to a lesser extent. The latent heat (Q_e) reached a maximum value during Gamma but was also very high during Delta, showing a large amount of seawater evaporation or vapor condensation (cloud formation) in the atmosphere for both events. Moreover, high values observed at the end of the time series (i.e., 25 October 2022) could be associated with the high air moisture caused by the floods left by the hurricanes. Ambient moisture condenses into raindrops, forming clouds and taking latent heat from the ocean's surface. Therefore, cloud formation could be the cause associated with the higher Q_e values observed at the end of this time series.

An empirical orthogonal function (EOF) analysis was performed on the alongshore component (u) of the ADCP currents, for the 3 years of the time series. Figure 4 shows this result from July to October 2020. Mode 1 accounted for 93 % of the explained variance; the time distribution depicted the strongest and fastest current fluctuation during the passage of these storms (orange lines in Fig. 4a.1), where negative (positive) values represent a westward (eastward) flow during Gamma (Delta). The spatial distribution (Fig. 4a.2) revealed a typical bottom Ekman layer distribution for the alongshore current, with a vertical shear $\left(\frac{\mathrm{\Delta }u}{\mathrm{\Delta }z}\right)$ of 0.005 s⁻¹ that is an alongshore current difference of 4 cm s⁻¹, in the 8 m spanned by the ADCP. During the passage of Gamma and Delta, $\frac{\mathrm{\Delta }u}{\mathrm{\Delta }z}$ presented a mean value of 0.024 s⁻¹ and a maximum of 0.077 s⁻¹. Figure 4b illustrates the alongshore current time–depth distribution. Coastal currents off Sisal flow preferentially toward the west, from the Caribbean Sea toward the Gulf of Mexico, 80 % of the time, considering the 3-year time series. However, it is common to find sporadic eastward flows, as shown by the red colors. During Gamma and Delta, the great mixing capabilities of these cyclones were evidenced by the strong and fast current fluctuation promoted in the water column over a few days (dashed lines in Fig. 4b).

Intense winds drove energetic waves (H_s>2 m) during the passage of both Gamma and Delta (Fig. 5a). Gamma induced NNW waves higher than 1.5 m and Tp>6 s from 3 to 6 October. Wave energy decreased during 7 October but increased by the end of the same day due to the passage of Delta, reaching H_s>2 m for a few hours due to the fast translation speed. The wave direction was from the NW (Fig. 5c), and the peak period was around 8 s (Fig. 5b). Therefore, the alongshore sediment transport is expected to occur in the eastward direction. After the passage of the tropical systems, the typical low-energy wave conditions associated with sea breezes were restored. The mean sea level at the coast increased 0.2 and 0.3 m for Gamma and Delta, respectively (Fig. 5d).

Figure 5 Time series of (a) significant wave height, (b) peak wave period, (c) mean wave direction measured at the ADCP deployed 10 km offshore, and (d) mean sea level (measured: solid line; predicted: dashed line).

4.2 Coastal impacts

4.2.1 Wrack and vegetation

Vegetation on the subaerial beach profile and foredune was present in Sisal before tropical storms Gamma and Delta (Figs. 6a and 7a). The more consolidated and dense vegetation was present 40 m away from the shoreline and remained unaltered after the meteorological events. However, high water levels associated with the passage of Gamma (2–5 October) either damaged or buried the beach vegetation in the low-elevation area in the vicinity of the port jetty and the central beach region (Figs. 6b and 7a). The further increase in the water levels and wave energy during the passage of Delta (8 October) affected all the pioneer vegetation on the beach's shoreward-most location (Fig. 6d).

Figure 6 Rectified images from the video monitoring system showing a (1) 1950 m and (2) 780 m stretch of Sisal beach for (a) 1 October 2020, (b) 2 October 2020, (c) 5 October 2020, and (d) 20 October 2020. Red boxes in (1) represent the area shown in (2). Images taken from http://tepeu.sisal.unam.mx/ (last access: November 2022).

Figure 7 (a) Wrack, (b) beach berm vegetation, and (c) flood distribution areas along the Sisal beach for 1, 2, 5, and 20 October 2020.

Prior to the storms, a small quantity of wrack was present along the beach (Figs. 6a and 7b). The increase in the incoming wave energy induced significant sediment transport and seabed erosion on the nearshore, causing the dislodgement of seagrasses. The seagrass was transported onshore by waves and currents and wrecked on the beach face, especially west of the pier (Figs. 6c and 7b). After the storms, the sea breezes redistributed the seagrass along the coast, accumulating most of the wrack east of the jetty (Figs. 6d and 7b).

4.2.2 Beach morphology and flooding

Beach profiles undertaken before and after Delta and Gamma allow assessing the storm impact on the beach morphology and estimating flooding. Pre- and post-storm beach profiles were analyzed to determine beach changes. The largest changes in both shoreline position and subaerial beach volume occurred in the vicinity of the pier and the port's jetty (see P02–P03 and P20–P22 in Fig. 8). The shoreline position east of the port (P20) and the pier (P02) retreated 22 and 10 m, respectively. On the other hand, 15 and 2 m shoreline advances were observed west of the jetty (P21) and the pier (P03). The mean shoreline change between transects P01 and P20 was −3 m, with transects P15–P18 showing a net increase. The latter suggests that significant eastward alongshore transport occurred during the storm sequence, induced by the NNW waves (Fig. 5c), redistributing the existing sediment east of the port (P19–P20) to adjacent transects (P15–P18). West of the port of Sisal (transects P21–P40) the mean shoreline advance was 4 m, with transects P27, P30, P31, P33, and P34 presenting a shoreline retreat in this area (Fig. 8a). The shoreline advance along P37–P40 seems to be related to the formation of a 0.20 m berm. It is important to point out that beach scarp erosion contributed to beach sediment accumulation at some transects.

Figure 8 (a, b) Shoreline and (c, d) subaerial beach volume changes obtained from the DGPS beach profiles east (P01–P20) and west (P21–P40) of the Sisal port.

The subaerial beach volume change, associated with cross-shore transport, presented an overall net increase (Fig. 8c and d). A lower impact on beach volume was observed at transects P04–P10 located in the area where significant seagrass wrack occurred. The sequence of storms did not induce a significant mean subaerial volume change (0.8 m³ m⁻¹), and the maximum volume increase (+11 m³ m⁻¹)/decrease (−18 m³ m⁻¹) corresponded to transects located west/east of the jetty (Fig. 8c and d). Significant volume losses also occurred at P02, P24, and P27. Away from the structures, significant sediment supply by the storms contributed to an increase in beach elevation (e.g., Tuck et al., 2021).

Beach morphology changes are also employed as a proxy for coastal flooding on the beach located in front of the coastal community of Sisal (i.e., P01–P20). Bed elevation increase was observed landward of the shoreline at most profiles (red and orange areas), whereas erosion was maximum east of the structures (blue areas) (Fig. 9a). The landward limit of observed bed changes was used as a proxy for the maximum horizontal swash excursion X_max (Fig. 9b), showing alongshore differences with a maximum closer to the port's jetty. The bed change associated with the maximum z was used as a proxy for the maximum water levels (Z_max = tide + storm surge + runup), implying that swash flows reaching that area were significant enough to induce sediment transport. The maximum elevation of such changes was found at P19 and the minimum at P07, corresponding to elevations of 1.7 and 0.5 m, respectively (Fig. 9c). The lower values are correlated with the areas that presented wrack coverage during the storms (Fig. 7b).

Figure 9 Coastal flooding derived from beach morphology changes east of the Sisal port.

The most vulnerable area to floods is located east of the jetty. The heavy rain of 2 October was capable of flooding a 200 m stretch of beach (Figs. 6b and 7c), and the succeeding forcings of Gamma propagated the flooded area farther to the east (Figs. 6c and 7c). The high vulnerability is due to dredging practices that have created a low-lying zone at this location.

4.2.3 Coastal aquifer

Monitoring wells W4, W5, and W7a located 20 km, 5 km, and 200 m from the coast provide information on the oceanic and terrestrial forcing on the coastal aquifer. Figure 10 shows the relative levels of the hydraulic head at each well and the precipitation at Sisal. Wells W7a and W5 show the diurnal tidal modulation on the hydraulic head (Fig. 10a). All wells show an increase in the water table owing to the recharge following the storms. This is more evident at the well located 20 km from Sisal (W4), which shows an increase of more than 2 m following the passage of Gamma and Delta and reaches a maximum level on 7 October. It is worth noting that coastal wells W5 and W7a do not show such an increase in the water table due to confined aquifer conditions that do not allow rapid infiltration of the precipitation (Fig. 10b); however, the storm effects can be seen in the change in amplitude (decrease) of the tide caused by the increase in aquifer discharge. The southern limit of such confinement is not well known, and hence back-barrier flooding might occur when the water table exceeds the confinement level south of the aquifer, preventing the hydraulic head in wells W7a and W5 from increasing further (Perry, 1989; Pino et al., 2011).

Figure 10 Times series of (a) the hydraulic head at three coastal monitoring wells (W4: solid black line; W5: solid gray line; W7a: red line), and (b) precipitation, recorded at the coastal weather station, in Sisal, Mexico. There are missing data due to sensor failure (a) during June 2021 and (b) from January to April and from September to October 2021.

4.3 Numerical modeling

4.3.1 Forecast modeling

Figures 11 and 12 show the time series of the measured and forecast data for 24, 48, and 72 h. Figure 11 shows that the 24 h forecasts adequately represent the variability in air temperature, wind direction, and atmospheric pressure. In the case of air temperature (Fig. 11a), the forecasts underestimate the diurnal variability, while during the storm events, they fit appropriately; thus, a relatively low correlation coefficient was obtained. For the wind speed, the forecasts consistently underestimate the magnitude of the wind (Fig. 11b), mainly during extreme events; however, the wind direction presents less bias than its magnitude (Fig. 11c). Regarding atmospheric pressure (Fig. 11d), the 24 h forecast shows a significantly high correlation coefficient, which suggests high reliability. Both the 48 and the 72 h forecasts fail to describe the atmospheric conditions during these events.

Figure 11 Measured (black line) and forecast modeled time series of (a) air temperature, (b) wind speed, (c) wind direction, (d) atmospheric pressure, and (e) precipitation. The red line indicates the 24 h forecast, the blue line the 48 h forecast, and the green line the 72 h forecast of WRF in Sisal, Yucatán, from 1 to 15 October 2020. The date format is mm/dd.

Figure 12 Measured (black line) and forecast modeled time series of (a) significant wave height, (b) wave direction, and (c) sea level anomaly. The red line indicates the 24 h forecast, the blue line the 48 h forecast, and the green line the 72 h forecast of WWIII for (a) and (b) and of the ADCIRC model for (c) in Sisal, Yucatán, from 1 to 15 October 2020. The date format is mm/dd.

In the same way as the atmospheric variables, the significant wave height, wave direction, and mean sea level were analyzed for the same period. Figure 12a shows that the forecasts adequately represent the temporal variability in the wave energy but significantly underestimate the significant height most of the time. The forecast that best fits the significant wave height variability is that of 24 h, mainly during Gamma. However, despite presenting a high correlation coefficient, the forecast significantly underpredicts the significant wave height during Delta. Regarding the wave direction, the three forecasts present an average bias of 1.6^∘ (Fig. 12b). Analyzing the change in sea level due to the approach of storms to the coast of Sisal, Fig. 12c shows that all the forecasts overpredict and underpredict the mean sea level during Gamma and Delta, respectively. The 24 h forecast is the one that best reproduces the sea level variability in the measured signal; however during extreme events, the 72 h forecast has a lower root-mean-square error (RMSE) and bias concerning the measured data than the 48 h forecast. Table 4 shows the statistical fit parameters that support the results obtained, and, additionally, Table 5 shows the correlation coefficient (CC), the RMSE, and the bias of each measured variable with respect to the forecasts.

Table 4 Statistical parameters at the 5 % significance level. df: degree of freedom; SD: standard deviation; CI: confidence interval.

Table 5 Representative statistics that validate the numerical simulations of the forecasts of the WRF, WWIII, and ADCIRC models. The correlation coefficient (CC), root-mean-square error (RMSE), and bias of the analyzed variables are shown. The standard deviation (SD) of the measured data is also shown.

4.3.2 Hazard assessment

The assessment of the waves and storm surges generated by Gamma and Delta is shown in Figs. 13–15. Regarding the waves generated by these events, Fig. 13 shows the results of the maximum wave height (Fig. 13a and c) and the maximum mean period (Fig. 13b and d) obtained during the entire path of Gamma (Fig. 13a and b) and Delta (Fig. 13c and d). As can be seen, the waves generated by the Delta event were much larger than those of Gamma (i.e., northern Yucatán Peninsula), indicating a greater hazard.

Figure 13 Envelopes of (a, c) maximum significant wave height and (b, d) maximum mean wave period, for both events: (a, b) Gamma and (c, d) Delta.

Figure 14 Maximum envelopes for the storm surge generated by events (a) Gamma and (b) Delta at Sisal, Yucatán.

Figure 15 (a) Subaerial volume change evolution for beach profiles located east of the port's jetty (storm period for Gamma and Delta: vertical dashed lines) and (b) time series of volume change for selected transects with different post-storm behavior (P7 – stable, P16 – accretion, P20 – erosion followed by recovery).

Figure 14 shows the hydrodynamic simulation results denoting the maximum attained storm surge levels for Sisal. Delta generated larger flood areas than Gamma, mainly affecting the settlements near the lagoon area. The hydrodynamic and wave models were initially run in forecast mode using the National Hurricane Center (NHC) prediction to forecast hazard areas and alert the local authorities. The results presented herein comprise the post-event analysis based on the setup used in the forecast model. The models used were calibrated for the Gulf of Mexico based on historical events but not for the Yucatán coast specifically. The storm surge assessment with the actual storm tracks shows higher values for Delta as it was a stronger event and passed closer to the study area. Sisal was barely affected by Gamma, while Delta created flooding in the deposition area updrift of the harbor and on the lagoon side of the town. In this sense, more studies need to be carried out in order to have a calibrated and validated forecast model for the area.

5 Discussions

5.1 Coastal resilience

Knowledge of coastal resilience to hurricane events is important for coastal communities on barrier islands. We analyze the time series of the beach morphology and the water table to determine how long the impact of these events lasted at the study site. Figure 15a shows the spatio-temporal evolution of the subaerial beach volume at transects located in front of the coastal community (i.e., P01 to P20) during the year following the storms. Observed morphological changes are significantly smaller at beach profiles located away from the structures (P05 to P10). The beach evolution is strongly influenced during the following months by the sediment impoundment by the port's jetty. Significant alongshore differences were observed in the subaerial beach evolution; hence, we present the time series at selected transects (Fig. 15b). For instance, the transect in the vicinity of the jetty (P20 in Fig. 15b) presented significant subaerial beach volume losses (20 m³ m⁻¹) after the passage of the hurricane events. The negative trend continues during the winter and reverses in the summer of 2021 without reaching the pre-disturbed condition. On the other hand, transect P16 shows an increase in the subaerial beach volume after the passage of the two events and remained relatively stable during the following year, while transect P7 did not present a significant subaerial volume change and remained in the pre-disturbed condition (Fig. 15b).

For the water table, the aquifer took about 6 months to reach the pre-disturbed conditions; this effect is more evident at well W7a. Previous studies have found similar results: the effects of hurricanes on the aquifer can last for months (Yam-Caamal and Graniel-Castro, 2014; Kovacs et al., 2017; Kovacs et al., 2017). Some of those effects are beneficial for the population because there is a greater recharge of the freshwater resource. Still, on the other hand, the flooding experienced inland (e.g., in the capital city) after hurricanes Gamma and Delta is caused partially by the storm drainage systems not being able to drain the stormwater efficiently (the water table is too high for the drainage to work properly); Canul-Macario (2022) estimated that this effect could last for more than 5 months.

5.2 Limitations

The present study focuses on obtaining high-resolution measurements that can help to calibrate numerical models to be further implemented in other areas. Nevertheless, the location of the hurricane's track can play an important role in the observed impact. Hence, concurrent measurements at different places along the northern Yucatán coast need to be considered in future field efforts. Field observations suggest that our DGPS measurement errors are significantly smaller than the observed changes in the beach and hence do not affect the conclusions reached in the present study regarding beach morphodynamics. On the other hand, the confinement of the aquifer (Perry, 1989; Villasuso Pino et al., 2011; Canul-Macario et al., 2020) plays an important role in the dynamics of the coastal aquifer. It is well known that the effects of the tide propagate further when compared to unconfined aquifers (White and Roberts, 1994; Canul-Macario et al., 2020). In the case of Yucatán, the confining layer boundaries are not well known; Villasuso Pino et al. (2011) show that its width decreases eastward. Therefore, more research is required on the coastal aquifer confinement to fully understand the coastal aquifer dynamics. Also, the effects of hurricanes on the aquifer of Yucatán are not well studied. Canul-Macario et al. (2020) studied the propagation of the atmospheric and meteorological tide on the coastal aquifer and found that meteorological tides propagate further inland. More detailed data are required to understand how hurricanes may impact the position of the saline interface in the aquifer, for example, and the role of the free aquifer on coastal flooding during such events (e.g., Geng et al., 2021), highlighting the need to estimate compound flooding at this location.

The numerical modeling presents limitations in implementing the forcing and boundary conditions due to the lack of some terrestrial and atmospheric processes not considered. The spatial and temporal resolution for the wind field associated with the storm passage is not high enough to capture the hydrodynamic response near the coast for such events that move/travel offshore. The uncertainty in the topography and bathymetry is important for the model implementation, since it controls wave transformation and coastal flooding. Moreover, field observations suggest that the flooding on the barrier island's lee side was associated with the high precipitation and the increase in the water table. These processes were not accounted for in the numerical modeling approach. Therefore, monitoring water levels in the wetlands would provide greater insights into the importance of precipitation in the flooding of the barrier island.

6 Conclusions

Tropical storms are important hazards on barrier islands of micro-tidal beaches along the northern Yucatán Peninsula. We investigate the impacts of hurricanes Gamma and Delta from the ocean to the coast using field observations and numerical models. Strong winds drive water mixing across an extensive area due to the shallow continental shelf. The study area is located more than ∼200 km away to the west of the center of the two storms. Although their impact lasted less than a week over the water column, the influence of both events in the oceanographic region was notable. Moreover, energetic waves and coastal currents induced significant sediment transport in the nearshore and were responsible for disaggregating macrophytes (seagrass) being further transported to the shore. The wrack's presence provided natural shoreline protection by increasing wave dissipation in shallow waters. Significant beach changes occurred due to the presence of coastal structures (e.g., port jetties), inducing alongshore sediment transport gradients. The beach located in the vicinity of the structure has a lower capability to recover from subaerial beach volume losses after the storms and did not reach the pre-storm condition after 1 year. Strong winds and low atmospheric pressure induced a storm surge on the order of the tidal range. The high water levels not only affected beach vegetation but also increased the subaerial beach volume due to cross-shore sediment transport. On the other hand, heavy rain increased the water level in the wetlands. The confinement of the aquifer plays an important role in the dynamics of the coastal aquifer and took several months to return to the pre-storm water level. More detailed data are required to understand how hurricanes may impact the position of the saline interface in the aquifer and the role of the free aquifer on coastal flooding during such events (e.g., Geng et al., 2021), highlighting the need to estimate compound flooding at this location.